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application  of  advanced  multivariate  inversion  techniques  to  generate  realistic,  comprehensive,  and  high-resolution 
3D  models  of  the  seismic  structure  of  the  crust  and  upper  mantle  that  satisfy  independent  geophysical  datasets.  Our 
focus  is  on  the  region  surrounding  Iran  from  the  east  coast  of  the  Mediterranean  in  the  west,  to  Pakistan  in  the  east, 
an  area  of  prime  importance  to  nuclear  explosion  monitoring  (NEM),  and  a  region  with  adequate  calibration  events 
to  validate  our  model  and  to  quantify  its  accuracy.  Specifically,  wc  plan  to  integrate  surface-wave  dispersion, 
receiver  function,  and  satellite  and  ground-based  gravity  observations  to  help  constrain  the  shallow  seismic  structure 
in  the  Arabian-Eurasian  collision  zone.  Building  on  our  earlier  work  combining  receiver  functions  and  surface  wave 
dispersion,  and  surface- wave  dispersion  and  gravity,  we  plan  to  continue  to  integrate  geophysical  data  sets  to  create 
more  compatible  earth  models.  We  also  intend  to  explore  geologically  based  smoothness  constraints  to  help  resolve 
sharp  features  in  the  underlying  shallow  3D  structure. 
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OBJECTIVES 

The  National  Nuclear  Security  Administration  (NNSA)  and  Air  Force  Research  Laboratory  (AFRL)  have  deeidcd  to 
investigate  3D  modeling  as  part  of  the  effort  to  improve  knowledge  of  Earth’s  congressional  and  shear  velocity 
structure  and  to  enable  a  reduction  of  uncertainty  in  our  ability  to  accurately  detect,  locate,  and  identify  small 
(mb  <  4)  seismic  events,  which  will  lead  to  improved  capabilities  for  NEM.  For  seismically  active  areas,  with  good 
ground-truth  event  coverage,  earth  models  with  limited  accuracy  can  be  corrected  by  interpolating  results  from 
nearby  ‘ground-truth’  events  (using  the  kriging  methodology)  making  it  possible  to  detect,  locate,  and  identify  large 
events  even  with  limited  resolution  of  Earth’s  structure.  However,  such  approaches  may  not  be  effective  for  smaller 
events,  and  remain  a  challenge  for  ascismic  regions.  To  improve  monitoring  capability  in  such  instances,  we  must 
develop  better  seismic  models. 

The  objective  of  this  work  is  to  help  improve  seismic  monitoring  technology  through  the  development  and 
application  of  advanced  multivariate  inversion  techniques  to  generate  realistic,  comprehensive,  and  high-rcsolution 
3D  models  of  the  seismic  structure  of  the  crust  and  upper  mantle  that  satisfy  independent  geophysieal  datasets.  Our 
foeus  is  on  the  region  surrounding  Iran  (Figure  1 )  from  the  cast  eoast  of  the  Mediterranean  in  the  west,  to  Pakistan  in 
the  cast,  an  area  of  prime  importance  to  NEM,  and  a  region  with  adequate  calibration  events  to  validate  our  models 
and  to  quantify  their  accuracy. 

Background 

Estimating  subsurface  geologic  variations  is  a  challenge.  Seismologists  have  worked  on  the  problem  for  more  than  a 
century  (c.g.,  Milne,  1899;  Maeelwane  and  Sohon,  1936;  Dahlen  and  Tromp,  1999).  As  data  quantity  and  quality 
and  computational  ability  have  improved,  we  have  made  important  advances  in  our  understanding  of  the  subsurfaee. 
Our  best  knowledge  applies  to  the  shallowest  regions  as  well  as  depths  with  the  sharpest  global  interfaces 
(sediment-basement  contacts,  the  base  of  the  crust,  base  of  the  mantle,  and  transitions  near  410  km  and  660  km 
depths),  where  resolution  is  improved  as  a  result  of  the  strong  interactions  of  body  w  aves  w  ith  sharp  geologic 
transitions  (e.g.,  Helmberger,  1968;  Langston,  1979;  Shearer,  1991;  Lay  t  al.,  2004).  We  have  also  done  well 
modeling  regions  with  smooth  velocity  changes  such  as  the  lower  mantle,  which  allows  us  to  exploit  the  information 
in  telcseismic  body- wave  travel  times  to  locate  seismic  sources  reasonably  well  (e  g.,  Kennctt  and  Engdahl,  1991). 
Still,  many  details  within  and  just  beneath  the  lithosphere  elude  us.  We  have  been  able  to  surmise  that  geologic 
variations  here  arc  substantial,  and  we  know  that  they  frustrate  attempts  to  use  robust  observations  such  as  regional 
seismic  travel  times  to  locate  events  in  many  parts  of  the  Earth  (e.g.,  Bondar  ct  al.,  2004). 

Travel-time  based  tomography  opened  the  doors  to  3D  imaging  but  the  models  remain  blurry,  often  suffer  from 
interpretational  ambiguity,  and  are  not  easily  used  to  predict  other,  independent  seismic  observations.  From  our  own 
analyses  (Maeeira  et  al.,  2005),  we  have  seen  how  high-resolution  surface-wave  tomography  fails  to  produce  the 
extremes  in  seismic  shear-wave  speed  that  are  evident  from  independent  observations  (we  diseuss  more  details 
below).  In  particular,  achieving  a  model  with  low  enough  seismie  wave  speeds  within  the  Tarim  Basin  to  match 
seismograms  from  high-quality  observations  remains  an  issue.  Waveform  tomography  methods  improve  the 
situation  somewhat,  including  information  from  both  the  amplitude  and  phase  of  the  signal,  but  restriction  of  these 
methods  to  lower  frequency  bands  limits  the  resolution  of  the  methods.  More  recent  finite-frequency  methods  (e.g., 
Zhou  et  al.,  2004;  Dahlen  and  Zhou,  2006)  and  adjoint  methods  (e.g.,  Tarantola,  1984;  Tromp  et  al.,  2008)  offer 
more  complete  approaches  to  modeling  waveforms  and  computing  sensitivity  kernels.  But  even  these  approaches 
face  limits  imposed  by  data  bandwidth.  In  any  event,  sueh  fully  3D  waveform  methods  could  benefit  greatly  from 
accurate,  if  approximate,  starting  models  derived  from  more  piecew  ise  interpretation  of  seismic  observables 
combined  with  other  observations. 

RESEARC  H  AC  COMPLISHED 

Our  award  was  made  late  in  the  Spring  of  2009,  so  much  of  this  section  is  a  review’  of  earlier  work  upon  which 
future  work  will  be  based  and  mueh  of  that  work  was  more  thoroughly  documented  in  earlier  Seismic  and 
Monitoring  Research  Review  Proceedings.  We  direct  the  reader  to  those  compilations  for  greater  detail.  We  begin 
writh  a  simple  conceptual  illustration  of  the  challenges  we  faee. 
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Including  Geologic  Information 

Geologic  boundaries  can  be  sharp,  and  a  sharp  geologic  transition  can  cause  problems  when  we  use  smooth  seismic 
models  to  predict  observations  for  events  or  stations  located  near  the  sharp  geologic  transition  (throughout  this 
report  sharp  should  be  interpreted  to  mean  sharp  relative  to  typieal  seismic  wavelengths).  A  smooth  model  can 
predict  the  travel  times  well  if  the  source  and  the  station  are  located  in  regions  with  smooth  velocity  structures  that 
arc  well  modeled  in  smooth  earth  models.  Thus  to  be  generally  applicable,  our  models  should  contain  sharp  features 
where  needed  (oeean-eontinent  transitions,  across  major  geologic  terrane  boundaries,  etc.).  How  can  we  reconstruct 
geologic  images  with  sharp  lateral  boundaries?  Or  how  can  we  use  long-period  observations  to  create  models  that 
we  can  use  to  estimate  short-period  wave  travel  times  and  amplitudes?  Imaging  sharp  boundaries  requires  broadband 
observations.  We  are,  however,  usually  limited  in  the  short  period  bands  beeause  of  their  strong  sensitivity  to 
structure  and  substantial  spatial  aliasing;  we  are  often  limited  at  the  long-period  end  by  noise  (i.e.,  long-period 
signals  from  many  small  events  are  smaller  than  background  seismic  ground  motions).  Even  when  they  are 
available,  long  period  signals  average  the  heterogeneity  over  broad  regions.  Thus  arises  the  conundrum  of  resolving 
geologic  detail  needed  to  explain  short-period  observations  using  long-period  data,  which  leads  to  the  problems 
associated  with  using  models  derived  using  long  periods  to  account  for  or  remove  propagation  effects  from 
short-pcnod  signals. 

To  produce  models  that  have  realistic  ‘sharp’  boundaries  requires  that  we  include  independent  information  on  the 
location  of  those  boundaries.  Such  information  is  available  (for  the  shallow  part  of  the  model)  in  independent  data 
sets  such  as  gravity,  surface  geologic  maps,  and  even  something  as  simple  as  topography.  As  part  of  this  work,  we 
plan  to  resolve  sharp  features  by  adapting  our  imaging  algorithms  to  allow  the  inclusion  of  geologie  information  on 
the  location  and  nature  of  the  boundary  into  shear-velocity  inversions  that  permit  such  features  (implemented 
through  custom  geologic  smoothness  constraints  that  allow  velocities  to  be  de-eorrelated  aeross  major  geologie 
transitions).  Including  a  priori  information  into  an  inversion  is  obviously  only  as  good  as  the  information  that  is 
included  Thus  the  inclusion  of  this  type  of  information  into  the  reconstruction  of  shear-veloeity  models  of  the 
subsurface  must  proceed  carefully  and  include  documentation  of  the  importance  of  the  assumed  a  priori  information 
on  the  resulting  model  From  a  probabilistic  approach,  uncertainties  on  the  a  priori  information  need  to  be  combined 
and  included  in  the  computations  of  the  shear-wave  speeds. 

Combining  Gravity  And  Rayleigh-Wave  Dispersion  Observations 

Inversion  of  surface  wave  dispersion  observations  is  a  standard  method  for  estimating  3D  shear  velocity  structure  of 
Earth's  crust  and  upper  mantle.  Nevertheless,  it  is  well  known  that  traditional  state-of-the-art  inversion  techniques 
suffer  from  poor  resolution  and  nonuniqueness,  especially  when  a  single  surface  wave  mode  is  used  (Huang  et  al., 
2003).  This  is  particularly  true  at  shallow  depths  where  the  shorter  periods,  which  are  primarily  sensitive  to  upper 
crustal  structures,  are  difficult  to  measure,  especially  in  tectonically  and  geologically  complex  areas.  On  the  other 
hand,  gravity  inversions  have  the  greatest  resolving  power  at  shallow'  depths  sinee  gravity  anomalies  decrease  in 
amplitude  and  increase  in  wavelength  with  increasing  depth.  Gravity  measurements  also  supply  constraints  on  rock 
density  variations.  Thus  by  combining  surface-wave  dispersion  and  gravity  observations  in  a  single  inversion,  we 
can  obtain  a  self-eonsistcnt  high-resolution  3D  shear-veloeity/density  model  with  increased  resolution  of  shallow 
geologic  structures. 

As  an  example,  consider  a  small  study  of  a  joint  surfaee-wave/gravity  inversion  performed  for  the  Tarim  Basin  in 
western  China,  which  shows  the  promise  of  this  approach  (Maeeira  and  Ammon,  2009).  We  used  gravity 
observations  extracted  from  the  global  gravity  model  derived  from  the  GRACE  satellite  mission  (Tapley  et  al., 

2005).  Specifically,  we  combined  Bouguer  gravity  anomalies  with  high-resolution  surface- wave  slowness 
tomographic  maps  (Maeeira  et  al.,  2005)  that  provide  group  velocity  dispersion  values  in  the  period  range  between 
8  and  100  s.  Figure  2  shows  the  gravity  (bottom  left)  and  dispersion  (top  left)  data  for  a  typical  cell  in  our  gridded 
model  together  with  the  fits  to  these  observations  resulting  from  the  inversion  of  only  dispersion  data  (second 
column  from  the  left),  from  the  inversion  of  only  gravity  data  (third  column  from  the  left),  and  from  the  joint 
inversion  of  both  data  (right  column).  The  best  fit  to  the  gravity  anomalies  comes  from  inverting  the  gravity  data 
alone,  meanwhile  that  same  model  is  not  able  to  fit  the  dispersion  observations.  In  the  same  way,  the  best  fit  to  the 
observed  dispersion  values  results  from  inverting  dispersion  observations  alone,  but  this  model  does  not  predict  the 
gravity  observations  adequately.  The  3D  velocity  model  obtained  from  the  joint  inversion  fits  both  data 
simultaneously  and  reasonably  well 

Figure  3  shows  the  3D  shear-wave  speeds  derived  from  the  joint  inversion  at  depths  of  2,  6,  10,  14,  28, 46,  55,  75, 
and  100  km.  In  general,  the  model  agrees  well  with  the  main  features  observed  in  previous  studies  (e.g.,  Villasenor 


5 


2009  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


et  al.,  2001 ;  Ritzwoller  et  al.,  2002;  Liang  ct  al.,  2004;  Hearn  et  al.,  2004;  Sun  et  al.,  2004;  Sun  and  Toksoz,  2006; 
Huang  and  Zhao,  2006).  At  shallow  depths,  low  velocities  in  the  Tarim  and  Junggar  basins  dominate  the  images. 
This  is  a  predictable  result  because  of  the  clear  presence  of  low  velocities  associated  with  the  known  major 
sedimentary  basins  in  the  tomographic  images  at  short  periods  (e.g.,  Maeeira  et  al.,  2005).  However,  the  seismic 
wave  speeds  in  the  model’s  basins  is  lower  than  those  obtained  without  including  the  gravity  in  the  inversion. 

We  think  that  this  is  an  improvement  of  our  3D  model  since  seismic  discrimination  has  proven  the  need  of  slower 
velocities  at  shorter  periods  (i.e.,  shallower  layers)  in  those  sedimentary  basins  (H.  J.  Patton,  see 
Acknowledgements).  Although  improved,  note  that  the  inversion  needs  further  refinement  since  the  high 
wavenumber  features  at  depths  of  75  and  100  km  must  be  artifacts  none  of  the  data  used  in  the  inversion  can 
constrain  such  features  (long-period  dispersion  and  gravity). 

To  quantify  the  improvement  in  the  shallow  parts  of  the  velocity  model,  we  tested  the  ability  of  the  3D  model  to 
predict  surfaec-wave  arrival  times  at  short  periods,  which  is  necessary  for  performing  surface-wave  magnitude 
measurements,  which  can  help  reduce  the  detection  threshold  for  seismic  discrimination  (Taylor  and  Patton,  2006). 
We  applied  the  method  described  by  Maeeira  (2006)  to  a  set  of  waveforms  from  26  nuclear  explosions.  The  new  3D 
model  is  able  to  better  predict  the  arrival  of  the  surface  waves  at  shorter  periods.  We  found  that  in  73%  of  the  eases, 
the  3D  model  from  the  joint  inversion  predicts  the  surface  wave  arrival  times  of  short  period  Rayleigh  waves  when 
the  dispersion-only  inversion  model  does  not.  The  combination  of  multiple  and  complementary  geophysical  data 
(surface  wave  dispersion  observations  and  gravity  measurements  in  this  ease)  in  a  simultaneous  inversion  not  only 
offers  a  simple  and  elegant  compromise  between  fitting  both  data  sets,  but  actually  improves  the  geophysical  model 
in  a  tangible  way  to  more  confidently  and  accurately  detect,  locate,  and  identify  small  seismic  events,  which  can 
help  improve  NEM  capabilities. 

CONCLUSIONS 


We  have  initiated  a  two-year  project  to  map  the  subsurface  geologic  variations  using  seismic  dispersion,  gravity,  and 
receiver- function  observations.  We  face  significant  challenges  in  our  efforts  to  include  effective  point  constraints  on 
structure  (receiver  functions)  with  the  spatially  continuous  surface-wave  tomography  and  gravity  observations.  Our 
work  complements  ongoing  work  at  Los  Alamos  National  Laboratory  to  integrate  body-wave  travel  times  into  the 
same  formalism.  The  basic  philosophy  is  that  models  that  explain  more  data  are  better.  The  ultimate  utility  of  the 
derived  earth  models  is  to  provide  improved  predictive  capabilities  for  routine  seismic  analyses  and  to  provide 
adequate  starting  models  for  3D  waveform  inversion  approaches. 
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Figure  1.  Map  of  focus  region  show  w  ith  topographic  and  bathymetric  shading  and  moderate  to  large 
earthquake  locations  (magnitudes  >  3.5  from  1990  to  Spring,  2008).  The  region  contains  the 
Arabian  plate  and  the  middle  segment  of  the  Alpine  to  Himalvan  collision  zone,  which  is 
constructed  primarily  of  Phanerozoic  terranes  amalgamated  onto  southern  Eurasia  during  the 
closing  of  the  Tethys  Ocean. 
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Figure  2.  Comparison  of  data  fits  from  the  inversion  of  only  surface  wave  dispersion  observations,  the 
inversion  of  only  gravity  observations,  and  the  joint  inversion  of  dispersion  and  gravity 
observations.  Top  panels  from  left  to  right:  surface  wave  dispersion  data  for  a  typical  cell  in  our 
gridded  model  (blue  line);  fit  to  the  dispersion  data  from  inverting  only  dispersion  data  (green 
line);  fit  to  the  dispersion  data  from  inverting  only  gravity  observations  (black  line);  fit  to  the 
dispersion  data  from  the  joint  inversion  (red  line).  Bottom  panels  from  left  to  right:  simple 
Bouguer  anomalies  for  the  region  under  study;  predicted  Bouguer  anomalies  from  the  model 
resulting  from  inverting  only  dispersion  observations;  predicted  Bouguer  anomalies  from  the 
model  resulting  from  inverting  only  gravity  observations;  predicted  Bouguer  anomalies  from  the 
model  resulting  from  the  joint  inversion. 
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Figure  3.  S-vvave  velocity  model  at  constant  depth  slices.  The  depth  of  each  image  is  shown  at  the  top  of 
each  map.  Velocity  values  arc  expressed  in  km/s.  Note  the  color  scheme  is  different  for  each 
image.  The  high  wavenumber  features  at  depth  arc  clearly  artifacts  of  the  simple  smoothing 
scheme  used  in  the  current  inversion  algorithm  -  none  of  the  data  (long-period  surface  wave 
dispersion  or  gravity)  can  constrain  such  sharp  features  at  those  depths. 
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